Optimization of plenoptic imaging systems

ABSTRACT

The spatial resolution of captured plenoptic images is enhanced. In one aspect, the plenoptic imaging process is modeled by a pupil image function (PIF), and a PIF inversion process is applied to the captured plenoptic image to produce a better resolution estimate of the object.

CROSS-REFERENCE TO RELATED APPLICATION(S)

This application is a continuation of U.S. patent application Ser. No.13/398,813, “Design and optimization of plenoptic imaging systems,”filed Feb. 16, 2012. The subject matter of all of the foregoing isincorporated herein by reference in its entirety.

BACKGROUND OF THE INVENTION

1. Field of the Invention

This invention relates generally to plenoptic imaging systems, and moreparticularly to reconstruction of images and image components inplenoptic imaging systems.

2. Description of the Related Art

The plenoptic imaging system has recently received increased attention.It can be used to recalculate a different focus point or point of viewof an object, based on digital processing of the captured plenopticimage. The plenoptic system also finds application in multi-modalimaging, using a multi-modal filter array in the plane of the pupilaperture. Each filter is imaged at the sensor, effectively producing amultiplexed image of the object for each imaging modality at the filterplane. Other applications for plenoptic imaging systems include varyingdepth of field imaging and high dynamic range imaging.

However, a plenoptic imaging system is not designed to capture a highfidelity image of the object. Rather, the “image” captured by aplenoptic system is intentionally scrambled in a manner that allows thecapture of additional lightfield information, but at the expense ofimage resolution. The loss of spatial resolution is a well-recognizedand debilitating drawback of plenoptic imaging systems. A lenslet arrayat the conventional image plane spreads light on to the sensor.Therefore each image of a given modality (such as each refocused image,or each wavelength image in multispectral imaging, etc.) is restrictedto as many pixels as there are lenslets.

Thus there is a need for approaches to combat this loss of spatialresolution.

SUMMARY OF THE INVENTION

The present invention overcomes various limitations by providingapproaches to enhance the spatial resolution based on captured plenopticimages. In one aspect, the plenoptic imaging process is modeled by apupil image function (PIF), which relates the object to the plenopticimage formed by that object. Given the PIF, a better resolution estimateof the object can then be calculated from the captured plenoptic imageby effectively inverting the PIF operation (although this may beachieved without actually inverting the PIF function).

In one implementation, the plenoptic imaging process is modeled as avector-matrix operation [Iv_(f)]=[PIFv][Iv_(o)], where [Iv_(o)] is avector representing the object (of different components in the object),[PIFv] is a matrix representing the PIF and [Iv_(f)] is a vectorrepresenting the captured plenoptic image. Often, the dimension of theobject vector Iv_(o) is much larger than that of the plenoptic imagevector Iv_(f), especially if the object vector Iv_(o) representsmultiple components within the object, for example from differentmodalities. That means the [PIFv] matrix is a rectangular matrix, forwhich no inverse matrix exists. However, the object vector can beestimated from the captured plenoptic image vector by solving a linearinverse problem, namely, computing an estimate [I*v_(o)], such that thereconstruction error ∥[PIFv][I*v_(o)]−[Iv_(f)]∥ is minimized. This willbe referred to as PIF inversion.

Such linear inverse problems can be formulated and solved in variousways. For example, different distance measures may be used for the errorcalculation. Constraints on the object vector may be included as part ofthe minimization process. For example, if the object vector is composedof different components of the object (e.g., red, green and bluecomponents of the object), the various components may have knownrelations to each other which could be enforced by using constraints inthe minimization process. Different optimization techniques may be usedto solve the minimization problem. Design choices of the matrix [PIFv]and selection of the optimization solving technique influence thequality of the estimated object vector [I*v_(o)]. In one approach, theindividual lenslets in the lenslet array in the plenoptic imaging systemare non-centrosymmetric in order to improve the behavior of theinversion problem.

The PIF can be estimated or approximated using different approaches. Awave propagation approach is more accurate and includes diffractioneffects. A geometrical optics approach is faster. The PIF may also beobtained by measurement of the optical response at the detector.

In another aspect, a resolution-enhanced plenoptic imaging systemincludes a conventional plenoptic imaging system coupled to digitalsignal processing that performs the PIF inversion calculation.

The basic PIF model described above can also be used for purposesbesides increasing the resolution of a captured plenoptic image. Forexample, the same basic principles can be used to assist in the designof these systems, to calibrate or test these systems, or to providereal-time automatic adjustments of these systems.

Other aspects of the invention include methods, devices and systemscorresponding to the concepts described above.

BRIEF DESCRIPTION OF THE DRAWINGS

The invention has other advantages and features which will be morereadily apparent from the following detailed description of theinvention and the appended claims, when taken in conjunction with theaccompanying drawings, in which:

FIG. 1 is a simplified diagram of a plenoptic imaging system.

FIGS. 2a-b illustrate the loss of resolution in a plenoptic imagingsystem.

FIGS. 3a-b illustrate spatial reconstruction in a plenoptic imagingsystem, according to the invention.

FIGS. 4a-e illustrate asymmetrical lenslets.

FIG. 5 illustrates portions of a PIFv matrix for an object with multiplespectral components.

FIG. 6 is a block diagram illustrating object reconstruction.

FIGS. 7a-d are block diagrams illustrating different modes usingfeedback based on error in the object estimate.

FIGS. 8a-b are block diagrams illustrating PIF inversion used forsubstance detection and for depth estimation, respectively.

FIGS. 9a-b illustrate determination of the PIF matrix using ray tracingthrough a filter module.

FIGS. 10a-b illustrate determination of the PIF matrix using ray tracingthrough a lenslet array.

The figures depict embodiments of the present invention for purposes ofillustration only. One skilled in the art will readily recognize fromthe following discussion that alternative embodiments of the structuresand methods illustrated herein may be employed without departing fromthe principles of the invention described herein.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

The figures and the following description relate to preferredembodiments by way of illustration only. It should be noted that fromthe following discussion, alternative embodiments of the structures andmethods disclosed herein will be readily recognized as viablealternatives that may be employed without departing from the principlesof what is claimed.

Sample Configuration.

FIG. 1 is a simplified diagram of a plenoptic imaging system. The systemincludes a primary imaging subsystem 110 (represented by a single lensin FIG. 1), a secondary imaging array 120 (represented by a lensletarray) and a sensor array 130. These form two overlapping imagingsubsystems, referred to as subsystem 1 and subsystem 2 in FIG. 1. Theplenoptic imaging system optionally may have a filter module 125positioned at a location conjugate to the sensor array 130. The filtermodule contains a number of spatially multiplexed filter cells, labeled1-3 in FIG. 1. For example, the filter cells could correspond todifferent modalities within the object.

The spatial coordinates (ξ, η) will be used at the object plane, (x, y)at the pupil plane for imaging subsystem 1, (u, v) at the plane of thelenslet array, and (t, w) at the sensor plane. The primary lens is adistance z₁ from the object, the lenslet array is a distance z₂ from theprimary lens, and the sensor is a distance z₃ from the lenslet array. InFIG. 1, the different components are each located at a single plane sothe distances between components is easy to define. However, in othersystems, the different components may be more complex (e.g., the primary“lens” may be a set of lenses), but the concept of the distances z₁, z₂,z₃ can be generalized using known principles. The focal lengths of theprimary lens and the lenslet array are represented by f₁ and f₂ andtheir diameters are D₁ and D₂, respectively. In addition, the differentcomponents do not have to be designed exactly as shown in FIG. 1. Forexample, the “primary lens” could be various combinations of elements,including lenses, mirrors and combinations of the two. Similarly, thesecondary imaging array could be a pinhole array, or a reflective array.

Ignoring the filter module 125 for the moment, in imaging subsystem 1,the object 150 is imaged by the primary lens 110 to produce an imagethat will be referred to as the “primary image.” This primary lens 110may be a camera imaging lens, microscope objective lens or any othersuch imaging system. The lenslet array 120 is placed approximately atthe location of the primary image. Each lenslet then images the pupil ofthe primary lens to the sensor plane. This is imaging subsystem 2, whichpartially overlaps with imaging subsystem 1. The image created at thesensor array 130 will be referred to as the “plenoptic image” in orderto avoid confusion with the “primary image.” The plenoptic image can bedivided into an array of subimages, corresponding to each of thelenslets. Note, however, that the subimages are images of the pupil ofimaging subsystem 1, and not of the object 150. In FIG. 1, the plenopticimage and subimages are labeled A1-C3. A1 generally corresponds toportion A of the object 150, as filtered by filter cell 1 in the filtermodule 125.

FIGS. 2a-b illustrate the loss of resolution in a plenoptic imagingsystem (without filter module 125). In this simulation, FIG. 2a showsthe object 150 and FIG. 2b shows the corresponding plenoptic image. Theconventional plenoptic image here has been obtained by simulating datacaptured behind every lenslet and binning all of that light into onepixel to give an in-focus image of the object. The plenoptic imagingsystem uses an 11×11 lenslet array 120. Note that the plenoptic imagebasically has an 11×11 resolution, determined by the lenslet array, eventhough there are many more than 11×11=121 sensors in the array 130. Thecaptured plenoptic image (at sensor 130) is not a high fidelityreconstruction of the object 150, even if the primary image incident onthe lenslet array 120 were a perfect reconstruction of the object. Theloss of spatial resolution is inherent to plenoptic imaging systems.

Image Reconstruction.

Even though this loss of spatial resolution is inherent to plenopticimaging systems, one of our insights is that the plenoptic image can beprocessed in a way to reconstruct a higher fidelity image of the object.For example, we will show that the plenoptic imaging system can bemodeled by

$\begin{matrix}{\begin{bmatrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}\begin{matrix}I_{f}^{1,1} \\I_{f}^{1,2}\end{matrix} \\\vdots\end{matrix} \\\vdots\end{matrix} \\I_{f}^{1,W}\end{matrix} \\I_{f}^{2,1}\end{matrix} \\\vdots\end{matrix} \\\vdots\end{matrix} \\I_{f}^{T,W}\end{bmatrix} = {\begin{bmatrix}{PIF}_{1,1}^{1,1} & {PIF}_{1,2}^{1,1} & \vdots & \vdots & {PIF}_{M,N}^{1,1} \\{PIF}_{1,1}^{1,2} & {PIF}_{1,2}^{1,2} & \vdots & \vdots & {PIF}_{M,N}^{1,2} \\\vdots & \vdots & \vdots & \vdots & \vdots \\\vdots & \vdots & \vdots & \vdots & \vdots \\{PIF}_{1,1}^{1,W} & {PIF}_{1,2}^{1,W} & \vdots & \vdots & {PIF}_{M,N}^{1,W} \\{PIF}_{1,1}^{2,1} & {PIF}_{1,2}^{2,1} & \vdots & \vdots & {PIF}_{M,N}^{2,1} \\\vdots & \vdots & \vdots & \vdots & \vdots \\\vdots & \vdots & \vdots & \vdots & \vdots \\{PIF}_{1,1}^{T,W} & {PIF}_{1,2}^{T,W} & \vdots & \vdots & {PIF}_{M,N}^{T,W}\end{bmatrix}\begin{bmatrix}I_{o,1,1} \\I_{o,1,2} \\\vdots \\\vdots \\I_{o,2,1} \\I_{o,2,2} \\\vdots \\\vdots \\I_{o,M,N}\end{bmatrix}}} & \left( {1A} \right)\end{matrix}$where I_(f) ^(T,W) is the intensity at sensor element T,W; I_(O,M,N) isthe intensity from object element M,N; and PIF_(M,N) ^(T,W) is theplenoptic imaging system response, which we refer to as the pupil imagefunction or PIF. T,W is the discretized version of sensor coordinates(t,w) and M,N is the discretized version of object coordinates (ξ,η).

Eq. 1A can be written more compactly as[Iv_(f)]=[PIFv][Iv_(o)]  (1B)where the additional “v” indicates that these are based on vectorizedimages. That is, the two-dimensional plenoptic image I_(f) isrepresented as a one-dimensional vector Iv_(f). Similarly, thetwo-dimensional object I_(o) is represented as a one-dimensional vectorIv_(o). Correspondingly, the four-dimensional PIF is represented as atwo-dimensional matrix PIFv.

The PIF can be calculated and approximated in different ways. We providetwo examples in the Appendix. Appendix A approximates a PIF based onwave propagation. Appendix B approximates a PIF based on geometricaloptics.

Based on Eq. 1B, an estimate [I*v_(o)] of the object can be calculatedby solving the following optimization problem:min_(Ivo)∥[PIFv][Iv _(o)]−[Iv _(f)]∥, subject to constraints on [Iv_(o)]  (2)The solution to this problem may be implemented as a closed-formsolution or via iterative algorithms, depending on the choice of thenorm (or combination of norms) ∥.∥ and the constraints on [Iv_(o)].

The following are some examples. For the problemmin_(Ivo)∥[PIFv][Iv_(o)]−[Iv_(f)]∥₂ ² with no additional constraints on[Iv_(o)], the solution is the least-squares solution[I*v_(o)]=[PIF*v][Iv_(f)], where [PIF*v] is the pseudo-inverse of[PIFv], and ∥.∥₂ is the Euclidean norm. For the problemmin_(Ivo)∥[PIFv][Iv_(o)]−[Iv_(f)]∥₂ ²+∥[Γ][Iv_(o)]∥₂ ² withregularization matrix Γ (Tikhonov regularization problem), the Tikhonovmatrix ([PIFv]^(T)[PIFv]−[Γ]^(T)[Γ])⁻¹[PIFv]^(T)Iv_(f) is the solution.For the problem min_(Ivo)∥[PIFv][Iv_(o)]−[Iv_(f)]∥₂ ² with upper andlower bounds on individual component of [Iv_(o)], a least square datafitting algorithm can be applied to calculate a solution. For thesparsity constraint regularization problemmin_(Ivo)∥[PIFv][Iv_(o)]−[Iv_(f)]∥₂ ²+λ∥[Iv_(o)]∥₀, solutions can beobtained via iterative optimization techniques such as Basis Pursuit orExpectation Maximization. For the L₁ regularization problem (LASSO)min_(Ivo)∥[PIFv][Iv_(o)]−[Iv_(f)]∥₂ ²+λ∥[Iv_(o)]∥₁ subject to∥[Iv_(o)]∥₁≦ε the solution can be obtained via convex optimizationtechniques. Other regularization techniques may be used.

Note that the observation vector [Iv_(f)] in Eq. (1B) is in general notthe result of a convolution of the system response matrix with theobject vector [Iv_(o)] and therefore direct deconvolution techniquessuch as Wiener filtering, etc. are not directly applicable.

The PIFv matrix may be obtained via simulation during system design, orby good calibration of each point in the object in a physical system, ormay be an approximation to this response. The PIFv matrix may also beupdated with techniques such as regularization, total variation, etc.,to incorporate approximations or errors in the knowledge of the PIFvmatrix.

Noise in the sensor may lead to the need for regularization in theinversion process. For a scientific sensor such as used in microscopeswe expect less noise. In one experiment, we obtained a reconstruction byusing Eq. 2 in an iterative least squares curve fit function in Matlabto estimate [Iv_(o)] from noisy sensor data [Iv_(f)]. We alsoincorporated a smoothing procedure in every iteration that removessharp, high frequency variations in the estimates of the objectreconstruction. Other variations of this procedure could be used.

FIGS. 3a-b illustrate spatial reconstruction, using the techniquedescribed above. These use the same object as FIGS. 2a-b . FIG. 3a showsreconstruction in a case without noise. FIG. 3B shows reconstructionusing iterative least squares, with SNR=40 dB.

In this example, the object is a 2D object. That is, the object distanceis the same for all points in the object. As a result, the whole objectshould be in-focus, and our corresponding reconstructions are alsoin-focus. However, the approach described above can also be used toobtain reconstruction for a 3D object for which the object plane isdefocused. In that case, some planes intersecting the object will bein-focus and some out-of-focus. If we consider a 3D object wheredifferent planes of the object may be in-focus or out-of-focus, then weupdate the PIFv matrix with the response of object points in x-y as wellas z depth. That is, the object vector [Iv_(o)] would contain multiplecomponents rather than a single object, each representing the componentof the object at a depth z. The PIFv matrix would also be expanded toinclude the contribution to [Iv_(f)] from each depth component. The PIFvmatrix then contains x-y-z calibration and its inversion would give datafor a 3D object. Depth information from other sources can be used toreduce the complexity of the inversion of this 3D response matrix.

Conditioning the PIFv Matrix.

Some techniques for calculating the solution of the PIFv inversionproblem may not always create useful results, even in the 2D, noise-lesscase. One common problem is that an erroneous twin-image solution candominate the reconstruction. The twin-image solution occurs when thereconstructed image is a combination of the object and its twin rotatedthrough 180 degrees.

This could be caused by symmetry in the response for some object points.In a radially symmetric optical system, a given object point almostalways has a diametrically opposite object point which gives a flippedversion of almost the same optical response (PIF). This means that thecoefficients in the PIFv matrix for one object point are almostidentical to the coefficients of another object point that is locateddiametrically opposite in the optical system at the plane of thelenslet. Thus, there is lack of enough diversity or randomness in thecoefficients in the PIFv matrix.

This can be addressed by incorporating asymmetry in the response for theobject points across a given lenslet. In one approach, we incorporateasymmetry in the shape/extent or amplitude or phase of the lenslets inthe lenslet array. FIGS. 4a-e illustrate some examples of asymmetricallenslets. In FIG. 4a , the shading represents the amplitude response ofthe lenslet. The absorption varies across the aperture of the lenslet.In FIGS. 4b-c , the shading represents the phase response of thelenslet. This asymmetric phase can be introduced by phase plates orother additional components, or could be the result of aberrationsdesigned into the lenslets. FIGS. 4d-e show two examples ofasymmetrically shaped lenslet apertures. In FIG. 4d , the aperture iscircular with a flat side. In FIG. 4e , the aperture is a regularpolygon with an odd number of sides. The lenslets have some variation(e.g., amplitude, phase or shape of aperture) that makes them lesssymmetric. All these cases are examples where the lenslet isnon-centrosymmetric (i.e., there is no 180 degree rotational symmetry).That is, the (n×n) pupil function matrix P for a lenslet satisfies thecondition P_(i, j)≠P_(n−i+1, n−j+1) for at least one pair (i, j) with1≦i, j≦n. Other forms of asymmetry could be a scratch or a mark orpattern across the face of a lenslet. This asymmetry imparts enoughvariation in the coefficients in the PIFv matrix to avert the twin-imagereconstruction error.

In one approach, the lenslets themselves are not altered. Rather, afilter or mask having amplitude or phase variations are placed in aconjugate object plane in order to implement such an asymmetry/variationin the response.

In a different approach, one could avoid using asymmetric lenslets andinstead use advanced image processing algorithms coupled to priorassumptions, such as object statistics or basis representations ofobjects, etc., to restrict the solution space for the PIF inversionproblem. However, this requires making assumptions about the object.These approaches can also be iterative and therefore computationallyexpensive.

Multiple Object Components

The description above can be used for reconstruction of a single object.In that case, Iv_(o) represents a single version of the object. Forexample, Iv_(o) might represent a grayscale image of the object, whereeach value in Iv, is the grayscale at a different pixel within theobject. The PIF matrix then represents the contribution of each pixelobject to the captured plenoptic image.

However, Iv_(o) and PIFv can also contain different components. Forexample, Iv_(o1) might be the component of the object at wavelength band1, Iv_(o2) might be the component of the object at wavelength band 2,and so on. The corresponding “component PIFs” could include PIFv₁, PIFv₂and so on, where the component PIFv₁ represents the contribution of thecomponent Iv_(o1) to the captured plenoptic image, component PIFv₂represents the contribution of the component Iv_(o2) to the capturedplenoptic image, etc. The basic equation [Iv_(f)]=[PIFv][Iv_(o)] stillholds, except that now the matrix [PIFv]=[[PIFv₁] . . . [PIFv_(K)]] andthe vector [Iv_(o)]=[[Iv_(o1)]^(T)[Iv_(o2)]^(T) . . .[Iv_(oK)]^(T)]^(T). Eqn. (1B) can then be written as[Iv_(f)]=[[PIFv₁][PIFv₂] . . . [PIFv_(K)]][[Iv_(o1)]^(T)[Iv_(o2)]^(T) .. . [Iv_(oK)]^(T)]^(T)  (3)The components could be based on wavelength, polarization, attenuation,object illumination or depth, for example. Using the notation of Eq. 1A,if Iv_(o) is a single object, then [PIFv] is a TW×MN matrix. If Iv_(o)is constructed from K components, then [PIFv] becomes a TW×KMN matrix.

In systems where a filter module 125 is used, the different componentsmight correspond to the different filter cells. In one exampleapplication, the filter cells have different spectral responses designedso that the filter module can be used to detect certain substances basedon these spectral responses. The object components Iv_(ok) are chosen asthe object contributions in each of the spectral regions. The PIFcomponents PIFv_(k) are then calculated as the system response for eachof the spectral regions. Given the captured plenoptic image Iv_(f), thePIF inversion then yields image estimates for each of the spectralcomponents Iv_(ok). The use of filter cells in the plenoptic imagingsystem facilitates this process since certain regions of the plenopticimage will be more strongly filtered by one spectral response thananother. Constraints on the object components may also be used duringthe PIF inversion. For example, there may be known correlation betweenthe different spectral components. Alternately, if there is anunfiltered object component, it might be equal to the sum of the otherspectral components. Note that Eq. 3 is general enough to account forspatial overlap between the different spectral components.

FIG. 5 illustrates components of a PIFv matrix. In this example, afilter module contains 3×3 filter cells, each with a different spectralresponse. The object Iv_(o) is divided into 9 corresponding spectralcomponents. The PIFv matrix then also has nine spectral components[PIFv_(k)], each of which is a TW×MN matrix. FIG. 5 shows each of thespectral components [PIFv_(k)], for a single object point (i.e. M=1,N=1) and for a simplified sensor with T=9, W=9. That is, the upper leftsquare shows the contribution to the plenoptic image (i.e., what thesensor detects) from the first spectral component of a single objectpoint. Most of the square is black because all but one of the filtercells blocks this spectral component. The next square shows thecontribution to the plenoptic image, from a second spectral component ofthe same object point. And so on, until the ninth square shows thecontribution from the ninth spectral component.

Other Modalities.

The model of Eqs. (1)-(3) can be used in a number of modes. Thedescription above described a situation where the plenoptic imagingsystem is used in an operational mode (as opposed to a calibration,testing or other mode). A plenoptic image is captured, and the goal isto reconstruct a high quality image of the original object (or of anobject component) from the captured plenoptic image. This will bereferred to as reconstruction mode. FIG. 6 is a block diagramillustrating a resolution-enhanced plenoptic imaging system used inreconstruction mode. An object 150 is incident on a plenoptic imagingsystem 510, which captures plenoptic image 520. The image captureprocess for the plenoptic imaging system is described by a pupil imagefunction (PIF) response. Signal processing 530 is used to invert thisprocess in order to obtain an estimate 550 of the original object/objectcomponent.

FIG. 7a is a block diagram based on the same PIF model, but used in adifferent manner. In this approach, the object 150 is known (orindependently estimated) and the estimated object 550 is compared to theactual object. An error metric 610 is calculated and used as feedback630 to modify the plenoptic imaging system 510 and/or the inversionprocess 530.

This general model can be used for different purposes. For example, itmay be used as part of the design process, as shown in FIG. 7b . Indesign mode, the plenoptic imaging system 510, captured plenopticimaging system 520 and signal processing 530 are computer models of thesystem being designed. The object 150 is also computer simulated. Theerror metric 610 is used to improve the design. It may also be used asfeedback 632 to improve the design of the plenoptic imaging system 510,the signal processing 530 or both. For example, positions, shapes andsizes of various elements in the plenoptic imaging system may beiterated based on the error metric. This includes the size and spacingof the lenslets and sampling at the sensor. The error metric may also beone component of a more inclusive metric. In another example, the filtermodule in the main lens may be translated in the plane perpendicular tothe optical axis such that the amount of light passing throughindividual filter cells is changes. In a third example, the position ofthe combined secondary imaging array and the sensor is changed along theoptical axis.

FIG. 7c shows an offline calibration mode. In this mode, the plenopticimaging system has been designed and built, but is being calibrated.This might occur in the factory or in the field. In this example, theobject 150 is a known calibration target. The error metric 610 is usedto calibrate 634 the already built plenoptic imaging system 510 and/orsignal processing 530. Example adjustments may include adjustments tothe physical position, spacing or orientation of components; to timing,gain, or other electronic attributes.

In FIG. 7d , the feedback loop is similar to FIG. 6c , except that itoccurs automatically in real-time 636. This would be the case forauto-adjustment features on the plenoptic imaging system.

In a last variation, FIGS. 7a-d are based on a metric that compares anestimated object 550 with the actual object 150. However, other metricsfor estimating properties of the object can also be used, as shown inFIGS. 8a-b . In addition, the PIF model can be used without having toexpressly calculate the estimated object.

In FIG. 8a , the task is to determine whether a specific substance ispresent based on analysis of different spectral components. For example,the PIF model may be based on these different spectral components, witha corresponding filter module used in the plenoptic imaging system.Conceptually, a PIF inversion process can be used to estimate eachspectral component, and these can then be further analyzed to determinewhether the substance is present. However, since the end goal issubstance detection, estimating the actual [Iv_(ok)] is an intermediatestep. In some cases, it may be possible to make the calculation 850 forsubstance detection without directly estimating the object components.The process shown in FIG. 8a can also be used in the various modalitiesshown in FIG. 7. For example, the system can be designed, calibratedand/or adjusted to reduce errors in substance detection (as opposed toerrors in object estimation). Errors can be measured by the rate offalse positives (system indicates that substance is present when it isnot) and the rate of false negatives (system indicates that substance isnot present when it is), for example.

Two examples of metrics that may be used for object classification arethe Fisher discriminant ratio and the Bhattacharyya distance. The Fisherdiscriminant ratio (FDR) can be used to evaluate the contrast betweentwo object classes. It is defined asFDR=( Y ₁ −Y ₂)^(T)(Σ_(Y1)+Σ_(Y2))⁻¹( Y ₁ −Y ₂)  (4)where Y ₁ and Y ₂ are the feature means of the two classes to bediscriminated and Σ_(Y1) and Σ_(Y2) are covariance matrices of the twoclasses.

The Bhattacharyya distance can be used to estimate classificationaccuracy based on a Bayes classifier. The Bhattacharyya distance is atheoretical value which is often used in an exponential function toestimate the upper bound of classification error. In one approach, ahigher order equation is used to estimate the classification error P_(e)based on the Bhattacharyya distance as shown below:P _(e) =a ₁ −a ₂ D _(B) +a ₃ D _(B) ² −a ₄ D _(B) ³ +a ₅ D _(B) ⁴ −a ₆ D_(B) ⁵  (5)where the coefficients are given by a₁=40.219, a₂=70.019, a₃=63.578,a₄=32.766, a₅=8.7172, and a₆=0.91875. The Bhattacharyya distance can becalculated using

$\begin{matrix}{D_{B} = {{\frac{1}{8}\left( {{\overset{\_}{Y}}_{1} - {\overset{\_}{Y}}_{2}} \right)^{T}\left( \frac{{\sum\limits_{Y\; 1}^{\;}{+ \;\sum\limits_{Y\; 2}^{\;}}}\;}{2} \right)^{- 1}\left( {{\overset{\_}{Y}}_{1} - {\overset{\_}{Y}}_{2}} \right)} + {\frac{1}{2}\ln\frac{\frac{{\sum\limits_{Y\; 1}^{\;}\;{+ \sum\limits_{Y\; 2}^{\;}}}\;}{2}}{{\sum\limits_{Y\; 1}^{\;}}^{1\text{/}2}{\sum\limits_{Y\; 2}^{\;}}^{1\text{/}2}}}}} & (6)\end{matrix}$

A support vector machine (SVM) can be used to estimate the hyperplaneseparating two classes that simultaneously minimizes the classificationerror and maximizes the geometric margin.

FIG. 8b gives another example, which is depth estimation. The task hereis to determine the depth to various objects. The object components arethe portions of the object which are at different depths. The PIFinversion process can then be used, either directly or indirectly, toestimate the different depth components and hence the depth estimates.This metric can also be used in different modalities. FIGS. 8a-b givejust two examples. Others will be apparent.

Although the detailed description contains many specifics, these shouldnot be construed as limiting the scope of the invention but merely asillustrating different examples and aspects of the invention. It shouldbe appreciated that the scope of the invention includes otherembodiments not discussed in detail above. Various other modifications,changes and variations which will be apparent to those skilled in theart may be made in the arrangement, operation and details of the methodand apparatus of the present invention disclosed herein withoutdeparting from the spirit and scope of the invention as defined in theappended claims. Therefore, the scope of the invention should bedetermined by the appended claims and their legal equivalents.

In the claims, reference to an element in the singular is not intendedto mean “one and only one” unless explicitly stated, but rather is meantto mean “one or more.” In addition, it is not necessary for a device ormethod to address every problem that is solvable by differentembodiments of the invention in order to be encompassed by the claims.

In alternate embodiments, portions of the invention can be implementedin computer hardware, firmware, software, combinations thereof and othertypes of processing power. Apparatus of the invention can be implementedin a computer program product tangibly embodied in a machine-readablestorage device for execution by a programmable processor; and methodsteps of the invention can be performed by a programmable processorexecuting a program of instructions to perform functions of theinvention by operating on input data and generating output. Theinvention can be implemented advantageously in one or more computerprograms that are executable on a programmable system including at leastone programmable processor coupled to receive data and instructionsfrom, and to transmit data and instructions to, a data storage system,at least one input device, and at least one output device. Each computerprogram can be implemented in a high-level procedural or object-orientedprogramming language, or in assembly or machine language if desired; andin any case, the language can be a compiled or interpreted language.Suitable processors include, by way of example, both general and specialpurpose microprocessors. Generally, a processor will receiveinstructions and data from a read-only memory and/or a random accessmemory. Generally, a computer will include one or more mass storagedevices for storing data files; such devices include magnetic disks,such as internal hard disks and removable disks; magneto-optical disks;and optical disks. Storage devices suitable for tangibly embodyingcomputer program instructions and data include all forms of non-volatilememory, including by way of example semiconductor memory devices, suchas EPROM, EEPROM, and flash memory devices; magnetic disks such asinternal hard disks and removable disks; magneto-optical disks; andCD-ROM disks. Any of the foregoing can be supplemented by, orincorporated in, ASICs (application-specific integrated circuits) andother forms of hardware.

In alternate embodiments, portions of the invention can be implementedby modules, which can take different physical forms. The term “module”is not meant to be limited to a specific physical form. Depending on thespecific application, modules can be implemented as hardware, firmware,software, and/or combinations of these. For example, the signalprocessing in FIGS. 6-7 could be implemented as dedicated circuitry(e.g., part of an ASIC), in order to take advantage of lower powerconsumption and higher speed. In other applications, the modules can beimplemented as software, typically running on digital signal processorsor even general-purpose processors. Various combinations can also beused. Furthermore, different modules can share common components or evenbe implemented by the same components. There may or may not be a clearboundary between different modules.

Depending on the form of the modules, the “coupling” between modules mayalso take different forms. Dedicated circuitry can be coupled to eachother by hardwiring or by accessing a common register or memorylocation, for example. Software “coupling” can occur by any number ofways to pass information between software components (or betweensoftware and hardware, if that is the case). The term “coupling” ismeant to include all of these and is not meant to be limited to ahardwired permanent connection between two components. In addition,there may be intervening elements. For example, when two elements aredescribed as being coupled to each other, this does not imply that theelements are directly coupled to each other nor does it preclude the useof other elements between the two.

APPENDIX A Example PIF Based on Wave Propagation

This derivation uses the nomenclature shown in FIG. 1 and thecorresponding text. Begin with imaging sub-system 1. Let the main lenshave a generalized pupil function denoted as P₁, with its correspondingimpulse response h₁:

$\begin{matrix}{{h_{1}\left( {u,{v;\xi},\eta} \right)} = {\frac{{\mathbb{e}}^{j\; k\; z_{1}}}{\lambda^{2}z_{1}z_{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{2}}\left( {u^{2} + v^{2}} \right)} \right\rbrack}{\exp\left\lbrack {\frac{j\; k}{2\; z_{1}}\left( {\xi^{2} + \eta^{2}} \right)} \right\rbrack}{\int{\int{{\mathbb{d}x}{\mathbb{d}y}\;{P_{1}\left( {x,y} \right)}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left( {x^{2} + y^{2}} \right)} \right\rbrack}\exp\left\{ {- {\frac{j\; k}{z_{2}}\left\lbrack {{\left( {u - {M\;\xi}} \right)x} + {\left( {v - {M\;\eta}} \right)y}} \right\rbrack}} \right\}}}}}} & ({A1})\end{matrix}$where λ is the wavelength of imaging,

$k = \frac{2\;\pi}{\lambda}$and the magnification from object to image plane is given by M=−z₂/z₁.Now consider the integral in Eq. (A1) as follows,

$\begin{matrix}{{\int{\int{{\mathbb{d}x}{\mathbb{d}y}\;{P_{1}\left( {x,y} \right)}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left( {x^{2} + y^{2}} \right)} \right\rbrack}\exp\left\{ {- {\frac{j\; k}{z_{2}}\left\lbrack {{\left( {u - {M\;\xi}} \right)x} + {\left( {v - {M\;\eta}} \right)y}} \right\rbrack}} \right\}}}} = {\int{\int{{\mathbb{d}x}{\mathbb{d}y}\;{P_{1}\left( {x,y} \right)}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left( {x^{2} + y^{2}} \right)} \right\rbrack}\exp\left\{ {{- j}\; 2\;{\pi\left\lbrack {{\left( {u - {M\;\xi}} \right)\frac{x}{\lambda\; z_{2}}} + {\left( {v - {M\;\eta}} \right)\frac{y}{\lambda\; z_{2}}}} \right\rbrack}} \right\}}}}} & ({A2})\end{matrix}$Substituting x′=x/λz₂ and y′=y/λz₂, the above term becomes

$\begin{matrix}{{{\lambda^{2}z_{1}z_{2}{\int{\int{{\mathbb{d}x^{\prime}}{\mathbb{d}y^{\prime}}{P_{1}\left( {{x^{\prime}\lambda\; z_{2}},{y^{\prime}\lambda\; z_{2}}} \right)}\exp\left\{ {\frac{j\; k}{2}{\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left\lbrack {\left( {x^{\prime}\lambda\; z_{2}} \right)^{2} + \left( {y^{\prime}\lambda\; z_{2}} \right)^{2}} \right\rbrack}} \right\}\exp\left\{ {{- j}\; 2\;{\pi\left\lbrack {{\left( {u - {M\;\xi}} \right)x^{\prime}} + {\left( {v - {M\;\eta}} \right)y^{\prime}}} \right\rbrack}} \right\}}}}} = {{{\lambda^{2}z_{1}z_{2}{{FT}\left( {{P_{1}\left( {{x^{\prime}\lambda\; z_{2}},{y^{\prime}\lambda\; z_{2}}} \right)}\exp\left\{ {\frac{j\; k}{2}{\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left\lbrack {\left( {x^{\prime}\lambda\; z_{2}} \right)^{2} + \left( {y^{\prime}\lambda\; z_{2}} \right)^{2}} \right\rbrack}} \right\}} \right)}}❘_{{{fx} = {({u - {M\;\xi}})}}{{fy} = {({v - {M\;\eta}})}}}} = {\lambda^{2}z_{1}z_{2}{h_{1}^{\prime}\left( {{u - {M\;\xi}},{v - {M\;\eta}}} \right)}}}},} & ({A3})\end{matrix}$where we have defined the term h′₁ for convenience as

$\begin{matrix}{{h_{1}^{\prime}\left( {u,v} \right)} = {{FT}\left( {{P_{1}\left( {{x\;\lambda\; z_{2}},{y\;\lambda\; z_{2}}} \right)}\exp\left\{ {\frac{j\; k}{2}{\left( {\frac{1}{z_{1}} + \frac{1}{z_{2}} - \frac{1}{f_{1}}} \right)\left\lbrack {\left( {x\;\lambda\; z_{2}} \right)^{2} + \left( {y\;\lambda\; z_{2}} \right)^{2}} \right\rbrack}} \right\}} \right)}} & \left( {A\; 4} \right)\end{matrix}$Eq. (A1) then reduces to

$\begin{matrix}{{{h_{1}\left( {u,{v;\xi},\eta} \right)} = {{\mathbb{e}}^{j\;{kz}_{1}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{2}}\left( {u^{2} + v^{2}} \right)} \right\rbrack}{\exp\left\lbrack {\frac{j\; k}{2\; z_{1}}\left( {\xi^{2} + \eta^{2}} \right)} \right\rbrack}}}{h_{1}^{\prime}\left( {{u - {M\;\xi}},\;{v - {M\;\eta}}} \right)}} & \left( {A\; 5} \right)\end{matrix}$

The field of the primary image formed at approximately the plane of thelenslet array isU _(i)(u, v)=∫∫dξdηh ₁(u,v;ξ,η)U _(o)(ξ,η)  (A6)where U₀(ξ,η) is the field of the object. Substituting Eq. (A5) yields

$\begin{matrix}{{{U_{i}\left( {u,v} \right)} = {{\mathbb{e}}^{j\;{kz}_{1}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{2}}\left( {u^{2} + v^{2}} \right)} \right\rbrack}{\int{\int{{\mathbb{d}\xi}{\mathbb{d}{{\eta U}_{0}\left( {\xi,\eta} \right)}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{1}}\left( {\xi^{2} + \eta^{2}} \right)} \right\rbrack}}}}}}{h_{1}^{\prime}\left( {{u - {M\;\xi}},\;{v - {M\;\eta}}} \right)}} & \left( {A\; 7} \right)\end{matrix}$Substituting ξ′=Mξ and η′=Mη in the above equation yields

$\begin{matrix}{{U_{i}\left( {u,v} \right)} = {\frac{{\mathbb{e}}^{j\;{kz}_{1}}}{M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{2}}\left( {u^{2} + v^{2}} \right)} \right\rbrack}{\int{\int{{\mathbb{d}\xi^{\prime}}{\mathbb{d}\eta^{\prime}}{U_{o}\left( {\frac{\xi^{\prime}}{M},\frac{\eta^{\prime}}{M}} \right)}\exp\left\{ {\frac{j\; k}{2\; z_{1}}\left\lbrack {\left( \frac{\xi^{\prime}}{M} \right)^{2} + \left( \frac{\eta^{\prime}}{M} \right)^{2}} \right\rbrack} \right\}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)}}}}}} & \left( {A\; 8} \right)\end{matrix}$

Now consider propagation from the lenslet array to the sensor. We assumeeach lenslet has a diameter D₂, focal length f₂, pupil function given byP₂ and there are M×N such lenslets in the array. The amplitudedistribution of the field after the lenslet array may be written as

$\begin{matrix}{{U_{i}^{\prime}\left( {u,v} \right)} = {{U_{i}\left( {u,v} \right)}{\sum\limits_{m}\;{\sum\limits_{n}\;{{P_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)}\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {u - {mD}_{2}} \right)^{2} + \left( {v - {nD}_{2}} \right)^{2}} \right\rbrack} \right\}}}}}} & \left( {A\; 9} \right)\end{matrix}$

Using Fresnel transformation to propagate this field to the sensor,which is located at a distance z₃ from the lenslet array, yields thefield of the plenoptic image U_(f)

$\begin{matrix}\begin{matrix}{{U_{f}\left( {t,w} \right)} = {\frac{{\mathbb{e}}^{j\;{kz}_{3}}}{j\;\lambda\; z_{3}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vU}_{i}\left( {u,v} \right)}}}}}}}}} \\{P_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)} \\{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {u - {mD}_{2}} \right)^{2} + \left( {v - {nD}_{2}} \right)^{2}} \right\rbrack} \right\}} \\{{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {u^{2} + v^{2}} \right)} \right\rbrack}{\exp\left\lbrack {\frac{{- j}\; k}{z_{3}}\left( {{ut} + {vw}} \right)} \right\rbrack}} \\{= {\frac{{\mathbb{e}}^{j\;{kz}_{3}}}{j\;\lambda\; z_{3}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;\exp}}}} \\{\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {mD}_{2} \right)^{2} + \left( {nD}_{2} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vU}_{1}\left( {u,v} \right)}}}}}} \\{P_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)} \\{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack} \\{\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}}\end{matrix} & \left( {A\; 10} \right)\end{matrix}$

Inserting Eq. (A8) into Eq. (A10) yields

$\begin{matrix}{{U_{f}\left( {t,w} \right)} = {\frac{{\mathbb{e}}^{j\;{kz}_{1}}{\mathbb{e}}^{j\;{kz}_{3}}}{{j\lambda}\; z_{3}M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {mD}_{2} \right)^{2} + \left( {nD}_{2} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vP}_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)}}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack}\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{\int{\int{{\mathbb{d}\xi^{\prime}}{\mathbb{d}\eta^{\prime}}{U_{o}\left( {\frac{\underset{\_}{\xi^{\prime}}}{M},\frac{\eta^{\prime}}{M}} \right)}\exp\left\{ {\frac{j\; k}{2\; z_{1}}\left\lbrack {\left( \frac{\underset{\_}{\xi^{\prime}}}{M} \right)^{2} + \left( \frac{\eta^{\prime}}{M} \right)^{2}} \right\rbrack} \right\}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)}}}}}}}}}}}} & \left( {A\; 11} \right)\end{matrix}$

Therefore the intensity of this field U_(f) at the sensor plane is givenby

$\begin{matrix}{{I_{f}\left( {t,w} \right)} = {\left\langle {{U_{f}\left( {t,{w;u},v,\xi^{\prime},\eta^{\prime}} \right)}{U_{f}^{*}\left( {t,{w;\overset{\sim}{u}},\overset{\sim}{v},\overset{\sim}{\xi^{\prime}},{\overset{\sim}{\eta}}^{\prime}} \right)}} \right\rangle = \left\langle {\frac{{\mathbb{e}}^{j\;{kz}_{1}}{\mathbb{e}}^{j\;{kz}_{3}}}{{j\lambda}\; z_{3}M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {mD}_{2} \right)^{2} + \left( {nD}_{2} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vP}_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)}}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack}\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{\int{\int{{\mathbb{d}\xi^{\prime}}{\mathbb{d}\eta^{\prime}}{U_{o}\left( {\frac{\xi^{\prime}}{M},\frac{\eta^{\prime}}{M}} \right)}\exp\left\{ {\frac{j\; k}{2\; z_{1}}\left\lbrack {\left( \frac{\xi^{\prime}}{M} \right)^{2} + \left( \frac{\eta^{\prime}}{M} \right)^{2}} \right\rbrack} \right\}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)} \times \left\langle {\frac{{\mathbb{e}}^{{- j}\;{kz}_{1}}{\mathbb{e}}^{{- j}\;{kz}_{3}}}{{- {j\lambda}}\; z_{3}M^{2}}{\exp\left\lbrack {{- \frac{j\; k}{2\; z_{3}}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{\overset{\sim}{m}}\;{\sum\limits_{\overset{\sim}{n}}\;{\exp\left\{ {\frac{{+ j}\; k}{2\; f_{2}}\left\lbrack {\left( {\overset{\sim}{m}D_{2}} \right)^{2} + \left( {\overset{\sim}{n}D_{2}} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}\overset{\sim}{u}}{\mathbb{d}\overset{\sim}{v}}{P_{2}^{*}\left( {{\overset{\sim}{u} - {\overset{\sim}{m}D_{2}}},{\overset{\sim}{v} - {\overset{\sim}{n}D_{2}}}} \right)}{\exp\left\lbrack {{- \frac{j\; k}{2}}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {{\overset{\sim}{u}}^{2} + {\overset{\sim}{v}}^{2}} \right)} \right\rbrack}\exp\left\{ {{+ j}\;{k\left\lbrack {{\overset{\sim}{u}\left( {\frac{t}{z_{3}} - \frac{\overset{\sim}{m}D_{2}}{f_{2}}} \right)} + {\overset{\sim}{v}\left( {\frac{w}{z_{3}} - \frac{\overset{\sim}{n}D_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{\int{\int{{\mathbb{d}{\overset{\sim}{\xi}}^{\prime}}{\mathbb{d}{\overset{\sim}{\eta}}^{\prime}}{U_{o}^{*}\left( {\frac{{\overset{\sim}{\xi}}^{\prime}}{M},\frac{{\overset{\sim}{\eta}}^{\prime}}{M}} \right)}\exp\left\{ {- {\frac{j\; k}{2\; z_{1}}\left\lbrack {\left( \frac{{\overset{\sim}{\xi}}^{\prime}}{M} \right)^{2} + \left( \frac{{\overset{\sim}{\eta}}^{\prime}}{M} \right)^{2}} \right\rbrack}} \right\}{h_{1}^{\prime*}\left( {{\overset{\sim}{u} - {\overset{\sim}{\xi}}^{\prime}},{\overset{\sim}{v} - {\overset{\sim}{\eta}}^{\prime}}} \right)}}}}}}}}}}} \right\rangle}}}}}}}}}} \right.}} & \left( {A\; 12} \right)\end{matrix}$

Considering an object creating spatial incoherence in the object field,

U _(o)(ξ′,η′)U* _(o)({tilde over (ξ)}′,{tilde over (η)}′)

=I_(o)(ξ′,η′)δ(ξ′−ξ′,η′−η′)  (A13)Substituting this into Eq. (A12) yields

$\begin{matrix}{{I_{f}\left( {t,w} \right)} = {{{\frac{{\mathbb{e}}^{j\;{kz}_{1}}{\mathbb{e}}^{j\;{kz}_{3}}}{{j\lambda}\; z_{3}M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}}}^{2}{\sum\limits_{m}\;{\sum\limits_{n}\;{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\left\lbrack {\left( {mD}_{2} \right)^{2} + \left( {nD}_{2} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vP}_{2}\left( {{u - {mD}_{2}},{v - {nD}_{2}}} \right)}}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack}\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\} \times {\sum\limits_{\overset{\sim}{m}}\;{\sum\limits_{\overset{\sim}{n}}\;{\exp\left\{ {\frac{{+ j}\; k}{2\; f_{2}}\left\lbrack {\left( {\overset{\sim}{m}D_{2}} \right)^{2} + \left( {\overset{\sim}{n}D_{2}} \right)^{2}} \right\rbrack} \right\}{\int{\int{{\mathbb{d}\overset{\sim}{u}}{\mathbb{d}\overset{\sim}{v}}{P_{2}^{*}\left( {{\overset{\sim}{u} - {\overset{\sim}{m}D_{2}}},{\overset{\sim}{v} - {\overset{\sim}{n}D_{2}}}} \right)}{\exp\left\lbrack {{- \frac{j\; k}{2}}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {{\overset{\sim}{u}}^{2} + {\overset{\sim}{v}}^{2}} \right)} \right\rbrack}\exp\left\{ {{+ j}\;{k\left\lbrack {{\overset{\sim}{u}\left( {\frac{t}{z_{3}} - \frac{\overset{\sim}{m}D_{2}}{f_{2}}} \right)} + {\overset{\sim}{v}\left( {\frac{w}{z_{3}} - \frac{\overset{\sim}{n}D_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{\int{\int{{\mathbb{d}\xi^{\prime}}{\mathbb{d}\eta^{\prime}}{I_{o}\left( {\frac{\xi^{\prime}}{M},\frac{\eta^{\prime}}{M}} \right)}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)}{h_{1}^{\prime*}\left( {{\overset{\sim}{u} - \xi^{\prime}},{\overset{\sim}{v} - \eta^{\prime}}} \right)}}}}}}}}}}}}}}}}}} & \left( {A\; 14} \right)\end{matrix}$This equation can be written as

$\begin{matrix}{{I_{f}\left( {t,w} \right)} = {\int{\int{{\mathbb{d}\xi^{\prime}}{\mathbb{d}{\overset{\sim}{\eta}}^{\prime}}{I_{o}\left( {\frac{\xi^{\prime}}{M},\frac{\eta^{\prime}}{M}} \right)}{\begin{matrix}{\frac{{\mathbb{e}}^{j\;{kz}_{1}}{\mathbb{e}}^{j\;{kz}_{3}}}{{j\lambda}\; z_{3}M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\begin{bmatrix}{\left( {mD}_{2} \right)^{2} +} \\\left( {nD}_{2} \right)^{2}\end{bmatrix}} \right\}}}}} \\{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vP}_{2}\begin{pmatrix}{{u - {mD}_{2}},} \\{v - {nD}_{2}}\end{pmatrix}}}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack}}}} \\{\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)}}\end{matrix}}^{2}}}}} & \left( {A\; 15} \right)\end{matrix}$

The assumptions related to incoherence allow us to separate the objectintensity I₀ and pull it outside the magnitude squared in the equationabove. The terms within the squared magnitude are all related toparameters of the optical system such as pupils, impulse responses forthe main lens and lenslet, etc. If this is considered the systemresponse, we can see that this term performs a separable operation withthe object intensity. Specifically, let

$\begin{matrix}{{PIF}_{\xi^{\prime},\eta^{\prime}}^{t,w} = {{{PIF}\left( {t,w,\xi^{\prime},\eta^{\prime}} \right)} = {\begin{matrix}{\frac{{\mathbb{e}}^{j\;{kz}_{1}}{\mathbb{e}}^{j\;{kz}_{3}}}{{j\lambda}\; z_{3}M^{2}}{\exp\left\lbrack {\frac{j\; k}{2\; z_{3}}\left( {t^{2} + w^{2}} \right)} \right\rbrack}{\sum\limits_{m}\;{\sum\limits_{n}\;{\exp\left\{ {\frac{{- j}\; k}{2\; f_{2}}\begin{bmatrix}{\left( {mD}_{2} \right)^{2} +} \\\left( {nD}_{2} \right)^{2}\end{bmatrix}} \right\}}}}} \\{\int{\int{{\mathbb{d}u}{\mathbb{d}{{vP}_{2}\begin{pmatrix}{{u - {mD}_{2}},} \\{v - {nD}_{2}}\end{pmatrix}}}{\exp\left\lbrack {\frac{j\; k}{2}\left( {\frac{1}{z_{2}} + \frac{1}{z_{3}} - \frac{1}{f_{2}}} \right)\left( {u^{2} + v^{2}} \right)} \right\rbrack}}}} \\{\exp\left\{ {{- j}\;{k\left\lbrack {{u\left( {\frac{t}{z_{3}} - \frac{{mD}_{2}}{f_{2}}} \right)} + {v\left( {\frac{w}{z_{3}} - \frac{{nD}_{2}}{f_{2}}} \right)}} \right\rbrack}} \right\}{h_{1}^{\prime}\left( {{u - \xi^{\prime}},{v - \eta^{\prime}}} \right)}}\end{matrix}}^{2}}} & \left( {A\; 16} \right)\end{matrix}$

This information about object separability and the lack of aconventional convolution with the system parameters iscounter-intuitive. However, following on this separability, the finalEq. (A15) for an incoherent plenoptic image can be re-written invector-matrix form to enable easy simulation. The term

$I_{o}\left( {\frac{\xi^{\prime}}{M},\frac{\eta^{\prime}}{M}} \right)$inside the integral refers to a single point in the object. The secondterm in the integral is represented by Eq. (A16). Thus, it can be seenthat this term is a function of the sensor plane coordinates (t, w) aswell as the single object point (ξ′, η′). When the entire integral isperformed over ξ′, η′ coordinates, we would effectively sum over all thepoints in an object and the result would be only a function of sensorcoordinates, (t, w), thus forming the basis for Eq. (1) of the maintext.

APPENDIX B Example PIF Based on Geometrical Optics

This derivation uses the nomenclature shown in FIG. 1 and thecorresponding text. For convenience, the following derivation isone-dimensional. The extension to two dimensions is straightforward. Thefollowing definitions are also used:

-   -   λ Design wavelength    -   D Diameter of primary lens 110    -   F Focal length of primary lens 110    -   d Diameter of microlens 120    -   f Focal length of microlens 120    -   k Lenslet counter    -   R Radius of curvature of concave microlens surface    -   t Thickness of microlens    -   n Refractive index of microlens    -   W Sensor size    -   Δ_(p) Size of sensor pixel    -   δ_(x) Sampling interval in the object plane    -   2M Number of object samples

One geometrics optics approach is based on ray trace matrices. The basicray trace matrix for the plenoptic imaging system of FIG. 1 is given by

$\begin{matrix}{{S\left\lbrack z_{3} \right\rbrack}{L\left\lbrack {- \frac{1}{f}} \right\rbrack}\left\{ t_{c_{k}} \right\}{S\left\lbrack z_{2} \right\rbrack}{L\left\lbrack {- \frac{1}{F}} \right\rbrack}\left\{ t_{p*} \right\}{S\left\lbrack z_{1} \right\rbrack}} & \left( {B\; 1} \right)\end{matrix}$where S[ ] is the matrix for free space propagation and L[ ] is thematrix for propagation through a lens. Referring to FIG. 1 and startingfrom the righthand side of expression B1, S[z₁] is the propagation fromthe object 150 to the primary lens 110. {t_(p*)} accounts for the lensaperture and filter module 125, if any.

$L\left\lbrack {- \frac{1}{F}} \right\rbrack$is propagation through the primary lens 110 and S[z₂] is propagationfrom the primary lens 110 to the lenslet array 120. {t_(c) _(k) }accounts for the aperture of the relevant lenslet in the array,

$L\left\lbrack {- \frac{1}{f}} \right\rbrack$is propagation through the lenslet and S[z₃] is propagation from thelenslet to the sensor. For convenience, wavelength dependence is notshown in Eq. B1. Note that there is some bookkeeping and coordinateshifting since some of the components are arrays of elements.Specifically, the filter module is an arrangement of filter cells andthe lenslet array is an array of lenslets. Thus, it must be determinedwhich cell in the filter module and which lenslet in the array istraversed by a ray before the correct ray matrices can be applied.

Eq. B1 is just one approach to tracing rays. Other ray-tracingapproaches can also be used. For example, exact ray tracing based onapplying Snell's law at each surface can also be used.

Regardless of the specific ray tracing method, ray tracing capabilitycan be used to estimate the PIF matrix. Recall that each value in thePIF matrix is the contribution of one object point/pixel to onepoint/pixel in the plenoptic image. That contribution can be estimatedby tracing rays.

One approach is based on Monte Carlo methods. A large number of rays istraced from an object point. The distribution of rays in angular spaceis representative of the energy distribution from the object. The numberof rays hitting different points in the plenoptic image is indicative ofthe contribution from that object point. This is repeated for manyobject points. In this way, the contribution from each object point toeach plenoptic image point (aka, the PIF matrix) can be estimated.

FIGS. 9a-b illustrate another approach based on identifying ray bundlesand which ray bundles hit different sensor elements. In this approach,each object point produces a ray bundle 910. These rays traverse throughthe system and eventually hit different sensor elements. By determiningwhat fraction of the original bundle 910 hits a particular sensorelement, we can then determine the corresponding elements in the PIFmatrix. This can be accomplished by tracing “border rays.” Border raysare the edge rays that define a ray bundle. The ray bundles are chosenso that the distribution of rays within the bundle can be determinedwithout expressly ray tracing them. The identification and tracing ofborder rays will sometimes involve forward ray tracing and sometimesbackward ray tracing.

Consider FIG. 9a . In this case, rays 911A and 911B are border rays forray bundle 910 since they define the edges of bundle 910. These edgerays are defined by the aperture of primary lens 110. However, when therays in bundle 910 propagate through filter module 125, not all the rayswill be treated the same since each filter cell 1, 2, 3 has a differentresponse. Accordingly, ray bundle 910 is subdivided into ray bundles920-1, 920-2 and 920-3, as shown by the dashed lines in FIG. 9b . Theseray bundles are defined by border rays 921C and 921D, in addition toborder rays 911A and 911B. All of these ray bundles 920 propagatethrough a single lenslet, so no further subdivision of the ray bundlesis required due to the lenslet array. Each ray bundle 920 hits thesensor array and the distribution of energy within each ray bundle 920can be calculated based on its coverage of the sensor array 130.

FIG. 10a shows an example where a ray bundle is subdivided because itpropagates through more than one lenslet. In FIG. 10a , there are threeray bundles 1020-1, 1020-2 and 1020-3 from the filter module. The solidlines are the border rays defining these three ray bundles. The middleray bundle 1020-2 hits the border of two lenslets and is split into tworay bundles 1020-2A and 1020-2B. The dashed lines are border raysdefining the extent of these two ray bundles.

FIG. 10b shows an example where a ray bundle intended for one lensletspills over into the subimage for an adjacent lenslet. In this example,lenslet 120A corresponds to sensor area 130A (denoted by the rectanglewith regions A1, A2 and A3), and lenslet 120B corresponds to sensor area130B. If there is clean separation between subimages, then all rays thathit lenslet 120A will be directed to sensor area 130A and contributeonly to the subimage produced by sensor area 130A. However, in thisexample, the rays in bundle 1020-3X actually spill over into sensor area130B. Thus, there is mixing of the subimages produced by sensor areas130A and 130B. This mixing is modeled in the PIF matrix and, therefore,can be unmixed by the PIF inversion process.

What is claimed is:
 1. A computer-implemented method comprising:accessing a design of a plenoptic imaging system; improving the designof the plenoptic imaging system by iteratively: modeling an operation ofthe plenoptic imaging system based on the design, the operation of theplenoptic imaging system producing plenoptic images of objects, theoperation of the plenoptic imaging system modeled by a pupil imagefunction (PIF) that relates an object to the plenoptic image formed bythat object; modeling reconstructions of the objects from the plenopticimages, said modeling based on applying a PIF inversion process to theplenoptic images; calculating a metric from the object reconstructions,the metric related to a spatial resolution of the objectreconstructions; and modifying the design of the plenoptic imagingsystem, such modification improving the calculated metric in a mannerthat improves the spatial resolution of the object reconstructions. 2.The method of claim 1 wherein the object reconstructions includedifferent depth components, and modifying the design comprises modifyingthe design to improve spatial resolution of the different depthcomponents.
 3. The method of claim 2 wherein the different depthcomponents include components located at different depth planes, andmodifying the design comprises modifying the design to improve spatialresolution at the different depth planes.
 4. The method of claim 3wherein modeling the operation of the plenoptic imaging system comprisesmodeling the operation for selected depth planes based on the design ofthe plenoptic imaging system, and modeling the operation for other depthplanes based on interpolating the models for the selected depth planes.5. The method of claim 1 wherein the object reconstructions includedifferent spectral components, and modifying the design comprisesmodifying the design to improve spatial resolution of the differentspectral components.
 6. The method of claim 1 wherein the objectreconstructions include different spectral components, and modifying thedesign comprises modifying the design to reduce mixing of the differentspectral components.
 7. The method of claim 1 wherein the objectreconstructions include different depth and spectral components, andmodifying the design comprises modifying the design to improve spatialresolution of the different depth and spectral components.
 8. The methodof claim 1 wherein the operation of the plenoptic imaging system ismodelled by a vector-matrix operation.
 9. The method of claim 1 whereinthe operation of the plenoptic imaging system is modelled based on wavepropagation applied to the design.
 10. The method of claim 1 wherein theoperation of the plenoptic imaging system is modelled based ongeometrical optics applied to the design.
 11. The method of claim 1wherein modifying the design is performed subject to object constraints.12. The method of claim 1 wherein modifying the design comprisesmodifying the design to reduce an error between the objectreconstructions and the objects.
 13. The method of claim 1 wherein thespatial resolution of the object reconstructions is higher than aspatial resolution of plenoptic images captured by the plenoptic imagingsystem.
 14. The method of claim 1, further comprising: constructing aplenoptic imaging system according to the modified design.
 15. Themethod of claim 1 wherein the object reconstructions include differentspectral components used to detect substances, and modifying the designcomprises modifying the design to improve the detection of thesubstances based on the different spectral components.
 16. Anon-transitory computer readable medium containing instructions that,when loaded into memory, cause a processor to carry out a methodcomprising: accessing a design of a plenoptic imaging system; improvingthe design of the plenoptic imaging system by iteratively: modeling anoperation of the plenoptic imaging system based on the design, theoperation of the plenoptic imaging system producing plenoptic images ofobjects, the operation of the plenoptic imaging system modeled by apupil image function (PIF) that relates an object to the plenoptic imageformed by that object; modeling reconstructions of the objects from theplenoptic images, said modeling based on applying a PIF inversionprocess to the plenoptic images; calculating a metric from the objectreconstructions, the metric related to a spatial resolution of theobject reconstructions; and modifying the design of the plenopticimaging system, such modification improving the calculated metric in amanner that improves the spatial resolution of the objectreconstructions.
 17. The non-transitory computer readable medium ofclaim 16 wherein the object reconstructions include different depthcomponents, and modifying the design comprises modifying the design toimprove spatial resolution of the different depth components.
 18. Thenon-transitory computer readable medium of claim 16 wherein the objectreconstructions include different spectral components, and modifying thedesign comprises modifying the design to improve spatial resolution ofthe different spectral components.
 19. The non-transitory computerreadable medium of claim 16 wherein the object reconstructions includedifferent depth and spectral components, and modifying the designcomprises modifying the design to improve spatial resolution of thedifferent depth and spectral components.
 20. The non-transitory computerreadable medium of claim 16 wherein the spatial resolution of the objectreconstructions is higher than a spatial resolution of plenoptic imagescaptured by the plenoptic imaging system.